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Abstract 

Angular momentum transport and particle acceleration during the magnetorotational instabil¬ 
ity (MRI) in a collisionless accretion disk are investigated using three-dimensional particle-in-cell 
(PIC) simulation. We show that the kinetic MRI can provide not only high energy particle accel¬ 
eration but also enhancement of angular momentum transport. We find that the plasma pressure 
anisotropy inside the channel flow with py > p± induced by active magnetic reconnection sup¬ 
presses the onset of subsequent reconnection, which in turn leads to high magnetic field saturation 
and enhancement of Maxwell stress tensor of angular momentum transport. Meanwhile, during 
the quiescent stage of reconnection the plasma isotropization progresses in the channel flow, and 
the anisotropic plasma with p± > py due to the dynamo action of MRI outside the channel flow 
contributes to rapid reconnection and strong particle acceleration. This efficient particle accelera¬ 
tion and enhanced angular momentum transport in a collisionless accretion disk may explain the 
origin of high energy particles observed around massive black holes. 
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A fundamental obstacle in our understanding of the universe is the need to explain the an¬ 
gular momentum redistribution in an accretion disk gravitationally rotating around a central 
object. It has been proposed that magnetorotational instability (MRI) is the most efficient 
mechanism for transporting angular momentum outward with inward mass motion BIZI- 
Based on magnetohydrodynamic (MHD) simulations [e.g.,|3H7j, it has been asserted that a 
weakly magnetized disk with an outwardly decreasing angular velocity gradient can provide 
angular momentum transport at a greatly enhanced rate by generating MHD turbulence. 

While the MHD framework is successful in explaining the “collisional” accretion disks, it 
is also important to study the dynamnics of “collisionless” accretion disks for some classes of 
astrophysical objects |H]. The accretion disk around the super-massive black hole Sagittarius 
A* at the center of our Galaxy is believed to be in a collisionless plasma state. This is because 
the accretion proceeds through a hot and low-density plasma in which the proton temper¬ 
ature is higher than the electron temperature [e.g., E]. In addition to the non-equilibrium 
temperature between protons and electrons, nonthermal high-energy particles are observed 

[e.g., [IUHI3]. 

Motivated by the observation of this collisionless accretion disk, Sharma et ah [mils] 
studied the MRI, including the effect of pressure anisotropy. Since the MRI involves the 
process of magnetic field amplification/dynamo, the perpendicular pressure is expected to 
be enhanced in the double adiabatic approximation |T6], which in turn modifies the MHD 
wave behavior. Meanwhile, because of the MHD waves generated by pressure anisotropy 
instabilities [e.g., E|, pressure isotropization occurs during MRI evolution. Riquelme et ah 
[I8] . Hoshino [I9] performed two-dimensional PIC simulations and confirmed the excitation 
of the mirror mode and the relaxation of pressure anisotropy studied by the previous fluid- 
based model [15]. In addition to the pressure anisotropy effect, the formation of a power-law 
energy spectrum during magnetic reconnection was pointed out. 

Although the previous two-dimensional PIC simulation showed the importance of the ki¬ 
netic accretion disk, the plasma transport process is generally different depending on whether 
it is considered in two-dimensional or three-dimensional space. The most important deter¬ 
minant of the efficiency of angular momentum transport in the collisionless MRI remains an 
open question. In this letter, we investigate for the first time the collisionless MRI using a 
three-dimensional PIC simulation, and argue that the angular momentum transport can be 
enhanced by the pressure anisotropy. 
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FIG. 1: Time evolution of the magnetorotational instability. Panels (a-b) show the magnetic field 
lines (greenish lines) and angular velocities in the background at Y = y/X = 1.91 and X = x/X = 
1.91 (color contour), and panels (c-e) depict the high density regions as reddish curved planes. The 
panels (b) and (c) are at the same time stage. Panel (f): The energy spectra during the MRI at 
^orbit = 0.31, 6.89, 7.18, 8.84 and 14.28. The dashed line is a Maxwellian fitting for Torut = 7.18. 

To study the kinetic accretion disk in three-dimensional space, we performed a PIC 
simulation in a local frame rotating with angular velocity VLqCz at a distance tq from the 
central object, and we include Coriolis, centrifugal, and gravitational forces in the equations 
of motion. We used tidal expansion of the effective potential with a constant q = —c^lnG/cllnr 
at ro, where g is 3/2 for a Keplerian disk. The equation of motion becomes 

du v 

— = e{E H— X B) — m7(2Go x v — 2qillxex)- (1) 

Our scheme was the same as that used in our previous two-dimensional MRI study, and we 
assumed that the local rotating velocity is much smaller than the speed of light [T^ . 
We adopted the shearing box boundary condition established by MHD simulations izi. 

For the initial condition, a drifting Maxwellian velocity distribution function was assumed 
in the local rotating frame with angular velocity f2o(^o)- The drift velocity in the ^/-direction 
Vy{x) was given by Vy{x) = rO(r) — rf2o(ro) — —qQo{rQ)x, and the radial velocity v^: and the 
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vertical velocity were both zero. In order to save CPU time, we set up the pair plasma, but 
the linear behavior of the MRI in the pair plasma was the same as that of ion-electron plasmas 
[T^ . A non-relativistic, isotropic plasma pressure with a high plasma (3 = 87i{p++P-)/Bq = 
1536 was assumed, where the electron and positron gas pressures were related to the thermal 
velocities Vt± by p± = {3/2)m±nv‘f^. The initial magnetic held was oriented purely vertical 
to the accretion disk; i.e., B = (0,0, Rq)- The ratio of the cyclotron frequency to the disk 
angular velocity was hxed at Uc±/^o = ±10, where = e±BQ/m±c. The grid size A was 
set to 2^^‘^(vt±/^p±), where /m± is the pair plasma frequency. The Alfven 

velocity is dehned as Va = B/y/87im±n, so that the plasma (3 is equal to 3v^^/V^. The 
parameters used were (Va/Uo)/A = 25, (nt±/Uc±)/A = 56.4, Va/c = 6.25 x 10“^. Aj,, 
Ny, and are the grid sizes in the x, y, and z directions, respectively, and we assumed 

Nx = Nz = Nz = 300 in this paper. = Ly = Lz = {NxA)/X = 1.91 is the physical size 

normalized by A = 2ttVa/^q- The number of particles per cell was set to Ap/cell = 40. 

Figure 1 shows the time evolution of the magnetic held lines (greenish lines) and the 
structure of the high density regions (sandwiched by the reddish curved planes). Color 
contours in the background at U = 1.91 and X = 1.91 show the angular velocity Vy in the 
local rotating frame. In the early stage at Torbit = ^ot/2TT = 0.31 in panel (a), the magnetic 
held lines are parallel to the z axis, and the Keplerian motion/diherential motion of Vy can 
be seen as the color contour at U = 1.91, where the reddish (bluish) region corresponds to 
a positive (negative) toroidal velocity. As time passes, the vertical magnetic helds start to 
get distorted due to the MRI, and they are stretched out in the toroidal direction because 
of the Keplerian motion at Torbit = 6.89 in panel (b). This stretching motion can amplify 
the magnetic held and form two inward- and outward-howing streams with a high plasma 
density and strong electric current, called the channel how. The reddish regions sandwiched 
by two surfaces in panel (c) show the high density channel how with p > (p) + 2ap where 
(p) and CTp are the average density and standard deviation of density distribution in the 
simulation domain, respectively. 

The amplihcation of the magnetic held stretched by the Keplerian motion may be bal¬ 
anced by the magnetic held dissipation caused by magnetic reconnection. Panel (d) at 
Torbit = 7.18 is the stage just after the onset of magnetic reconnection, and the break of 
the laminar channel how seen at Torbit = 6.89 can be observed. After the hrst onset of 
reconnection, subsequent reconnection occurs intermittently in several diherent sites in the 
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turbulent channel flows, and the formation of the channel flow with a strong magnetic held 
by MRI dynamo and destruction by reconnection occurs repeatedly. 

Panel (f) shows the time evolution of the energy spectra, where the horizontal and vertical 
axes are the particle energy normalized by the rest mass energy and number density N{e). 
Before the first onset of magnetic reconnection at Torbit = 6.89, the plasmas are gradually 
heated from the initial cold Maxwellian plasma. After the onset of reconnection at Torbit = 
7.18, we can clearly observe nonthermal particles above e/mc^ > 0.2. The dashed line, for 
reference, is a Maxwellian spectrum fitted by T/mc^ = 0.121. The nonthermal population 
continues to grow, and the spectrum can be approximated by a single power law function 
with N{e) oc at Torbit = 8.84. After Torbit = 9 ~ 10, the spectrum slop becomes 

harder in the high energy range from e/mc^ ~ 10 to 10^. The spectrum hardening might 
be due to the stochastic, multiple reconnection process IZH, but note that the maximum 
attainable energy in the system, whose gyro-radius is almost the same as the simulation box 
size, is e/mc^ ~ 10^. Then the spectrum deformation might be related to the accumulation 
of high energy particles around the maximum attainable energy. As already discussed by 
the previous 2D PIC simulations [HI HE], the pressure anisotropy with > p\\ is generated 
in our 3D simulation by the MRI dynamo (see Figure 3), which can contribute to rapid 
reconnection [22] and particle acceleration |19j . 

Let us take a look at the history of kinetic and magnetic held energies in the top panel of 
Figure 2. The energies are normalized by the initial magnetic held energy. As time goes on, 
both the kinetic and magnetic held energies increase, but the rapid increase of the magnetic 
held energy can be observed at around Torbit ~ 6, and the instantaneous plasma (3 becomes 
of the order of unity. Our PIC simulation in the local rotating system has been carried 
out using the open shearing box boundary condition [20] , and the plasmas can acquire their 
energies by accretion toward the center of gravity. Around Torbit ~ 8, the total magnetic held 
energy reaches to its maximum, and then starts to decrease until Torbit ~ 9- After Torbit ~ 9, 
both the magnetic held and kinetic energies remain almost constant with fluctuations. 

This time evolution is basically similar to our previous 2D PIC simulation result [19] . In 
the early nonlinear stage, we observe the formation of channel flows with the amplification 
of the magnetic held and the subsequent break of the channel flow by reconnection. The 
main difference, however, is seen in the late nonlinear stage with the turbulent/intermittent 
reconnection after Torbit ^ 9. In our previous 2D simulation, after the first onset of recon- 
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FIG. 2: (Top) Time evolution of the kinetic and magnetic field energies. (Bottom) The a parameter 
(black) and the contribution of Reynolds stress wr/p (green), Maxwell stress wm/p (red), and 
anisotropy stress waIp (blue). 

nection, the channel flows are destroyed, and a conple of large magnetic islands are formed 
in the system. In this 3D simnlation, however, the channel flows are preserved beyond the 
first onset of reconnection, and magnetic reconnection occurs in several different locations 
with the dynamic motion. The turbulent/intermittent reconnection in the late phase can 
be seen in 3D MHD simulations as well |3H7]. 

The most intriguing result in our kinetic MRI simulation is the enhancement of the 
angular momentum transport. The bottom panel in Figure 2 shows the time evolution of 
parameter a, which is used in the standard accretion disk model [23] and can be defined as 
« = Wxy/p, where p and w^y are the volume-averaged, instantaneous plasma pressure, and 
stress tensor, respectively. The stress tensor Wxy, which is related to the energy dissipation 
rate in the system, can be calculated as follows: 

Wxy = pVx{vy -b gDox) - -b (2) 

The terms on the right-hand side represent the Reynolds (wr), Maxwell {wm), and 
anisotropy {wa) stresses, respectively |15|. During the active reconnection phase between 
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7 < Torbit < 9, we found that a reached 0(1) with wm > wr ~ wa-, and during the late 
stage of Torbit > 9, a 0(10 with vjm ~ Wa > wr, which suggests a much more efficient 
angular momentum transport than the one discussed previously with a 0(10-3 - 10-2) 
[3 ED]. Note that a in some simulation results were normalized by the initial pressure po 
instead of the instantaneous pressure p. 

Since the parameter a is approximated by, 



(3) 


the enhancement of parameter a is related to a higher saturation of the magnetic field, 
which can be determined from the balance between the magnetic held amplihcation due to 
the MRI dynamo and magnetic held dissipation by reconnection. The dissipated magnetic 
held energy is deposited as thermal energy. If the onset of reconnection requires a high 
magnetic held in the collisionless system, then the plasma (3 becomes small and a large a 
can be realized. 

To understand the dynamics of reconnection, we focus on the pressure anisotropy. Panels 
(a) and (b) in Figure 3 at Torbit = 14.28 show, respectively, two dimensional (x, z) slices 
of plasma density and pressure anisotropy pu/py at the position Y = 0.96. The higher 
density regions around Z = 0.64 and 1.60 correspond to the so-called channel hows where 
the magnetic held polarity changes. At the center of the channel how, one can see p±/p|| < 1, 
while p±/p|| > 1 for the other regions. Panel (c) shows the relationship between the plasma 
density and p±/p||, and panel (d) is the histogram of occurrence frequency. This behavior can 
be basically understood using the double adiabatic equation of state with p_\_/pB = const. 
and p||T2/p3 _ ^onsf. [151 [IH ESj- The production of p±/p\\ > 1 is simply due to the 
magnetic held amplihcation of the MRI dynamo, and the formation of p\\/p± > 1 is because 
of magnetic reconnection at the center of the channel how, where the total B is dissipated 
while the plasma density is compressed. In the kinetic perspective, it is known that the 
pressure anisotropy can be produced by the Alfvenic beams along the plasma sheet boundary, 
which are emanating from the magnetic dihusion region m- 

The onset of magnetic reconnection is still a controversial issue, but the linear growth 
rate of the collisionless tearing mode under the pressure anisotropy would be sufficient for 
our argument. This is given, for example by simplifying Eq.(40) in [22]; as follows: 
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FIG. 3: (a) A slice of the three-dimensional density in x — z plane at y = 0.96, and (b) a slice 
of the pressure anisotropy p±/p\\ in the same plane. The color scales are shown on the right, (c) 
The color contour of the occurrence frequency in {p,p±/p\\) plane, and (d) the histogram of the 
occurrence frequency as a function of p±/p\\ in the logarithmic scale. 


where k, S, Vg, and Vth are the wave number, the thickness of the current sheet, gyro-radius, 
and thermal velocity, respectively. At the saturation stage, Vth < c, B/Bq ~ 230, and 
5/A > 10. We then obtain the estimates of < 0.089 and (1 — k‘^5‘^)/k5 ~ 0(1). 

On the other hand, the pressure anisotropy is 1 —p±/p\\ < 0.5 from Figure 3 (d). Therefore, 
it is highly possible that the successive reconnection in the channel flow is suppressed by 
anisotropic plasma of p±/p\\ < 1, which is formed by the preceding reconnection. The kinetic 
magnetic reconnection involves a deterrent effect to the successive magnetic dissipation, and 
as a result, the high magnetic held is realized before the onset of reconnection. 

To conhrm the enhancement of a parameter by the p\\ > p±_ effect in the channel how, we 
performed another simulation including an isotropization model for an anisotropic plasma 
after the formation of the channel how, Torbit > 6.89. In this model, we added a weak 
external random magnetic held 5B in Eq.(l) only for the channel how region which is 
roughly characterized by the weak magnetic held region of \B/Bq\ < 50. Namely, we used 
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FIG. 4: (Top) Time evolution of the volume integrated pressure anisotropy. The dashed and solid 
lines show the cases with and without the isotropization model, respectively. (Bottom) The a 
parameter (black) and the contribution of Reynolds stress wr/p (green), Maxwell stress wm/p 
(red), and anisotropy stress waIp (blue). 


the equation of motion given by. 


^ = e (e +- X {B + SB)] + (other forces). 


(5) 


c ' / 

where we assumed the white noise SB with |(5i?|/|i?| = 2.5. 

Shown in the top panel of Figure 4 is the time history of the volume integrated pressure 
anisotropy p±/p\\- The dashed and solid lines correspond to the cases with and without the 
above isotropization model, respectively. As we expected, one can hnd that the p±/p\\ with 
the isotropization model in the channel flow is larger than that withont the isotropization, 
because the anisotropic plasma with py > p± in the channel flow can be reduced by the 
isotropization model. 

Under this weak isotropization in the channel flow, let us study the time history of the 
a parameter and the Reynolds wr/p, Maxwell wm/p, and anisotropy stress wa/p with and 
withont the isotropization model in the bottom panel. The dashed and solid lines show 
the cases with and without the isotropization model, respectively. One can hnd that a 
parameter can be rednced nnder the isotropization model. The magnitnde of the anisotropy 
stress wa/p does not change between two cases with and without the isotropization model, 
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while the change of the Maxwell stress wm/p becomes large under the isotropization model. 
This suggests that the isotropization in the channel flow plays an important role on the 
magnetic field generation during the MRI dynamo. 

In summary, we have investigated for the first time a three-dimensional, collisionless 
MRI in a local rotating system, and have shown that an anisotropic pressure of p\\/p± > 1 is 
maintained in the channel flow during the MRI, which leads to high magnetic field saturation 
and an enhanced a parameter. During the quiescent stage of reconnection the isotropization 
of the anisotropic plasma progresses in the channel flow. After the plasma isotropization 
the anisotropic plasma with p±/p\\ > 1 outside the channel flow region may contribute to a 
rapid reconnection and nonthermal particle generation [T9] . 
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